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ABSTRACT 


A transient overpower (TOP) accident in a Liquid Metal 
Fast Breeder Reactor (EMFBR) is considered. The analysis 
is formulated to modei the dynamic response of the reactor 
fuel subassembly during the initial period of the postulated 
Overrower transient. An equivalent cylindrical cell is used 
to model the fuel subassembly, The governing neutronic and 
heat transpert equations for each region (fuel, clad, and 
coolant) of the equivalent cylindrical cell are developed. 
Nuclear Doppler broadening feedback is included in the dynam- 
ic model making the coupled equations non-linear. The re- 
suiting non-linear partial differential field equations are 
transformed into a system of ordinary differential equations 
by the finite element method. An isoparametric, quadratic, 
rectangular element is used for the discretization of the 
spatial domain. When using the finite element method, large 
system matrices may result. To facilitate selution of these 
large systems, an optimum compacting scheme is utilized. 
The implicit Gear's method is used for the solution of the 
system of ordinary differential equations. The results for 
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I. INTRODUCTION 


As the world's fossi? fuel resources are depleted, more 
emphasis {is being placed on the breeder reactor as a poten- 
tial means of solving the coming energy crisis. While the 
development of new energy sources is being pushed, equal ef- 
fort is heing given to the maintenance of an environmental- 
ly clean world. To this end, the safety of breeder reactors 
fs receiving a considerable amount of attention before assum- 
ing that the breeder reactor is the answer to the energy 
problem, : 

The Liquid Metal Fast Breeder Reactor (LMFBR) appears to 
be one cf the most promising breeder reactors. Most engineers 
will concede there is little probability of a suclear explo- 
sion occurring in the operation of a nuclear reactor. If 
major concern to engineers is the Joss of coolant accident 
(LOCA) and the transient overpower accident (TOP). The pres- 
ent analysis is concernec with a TOP accident in a LMFBR. 

The analysis is formulated to model the dynamic response of 
the reactor fue] subassembly during the initial period of the 
postulated overpower transient. The primary consideration is 
given to the early response of this fuel subassembly to 
various conditions of disturbances. The phenomenon which oc- 
curs after core disassembly (i.e., clad melting) is not the 
concern of this analysis. Only the time prior to clad melting 


is being considered, 
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No consideration is given here as to how the overpower 
transient occurs or to why the safety features of the reac- 
tor did not operate properly. It is postulated that the 
accident has occurred. In this analysis, the TOP accident 
1s created by either a step increase in reactivity, a ramp 
increase in reactivity, or a combination of both. 

An inherent safety feature of most reactors, nuclear 
Doppler broadening feedback, is included in the dynamic model 
of the fuel subassembly. The Doppler feedback acts to reduce 
the effect of the excursion. Consideration of this feedback 
creates @ non-linear system model which is described by a 
non-linear, initial-boundary-value problem. 

The conventional method of solution uses the standard 
point kinetics formulation. Recent studies have pointed out 
a non-negligible error in this model [1], particularly with 
asymmetric disturbances [2], or space-dependent feedback [3]. 
In Ref. [4], a somewhat novel approach of using the finite 
element method (FEM) for the space-time dependent solution 
of the reactor dynamics problem was demonstrated. The FEM 
is effective in handling these asymmetric disturbances and 
Space-dependent feedbacks. Therefore, the finite element 
method was used sc that the spatial effects on the postu- 
Tated problem may be studied turther. 

The purpose of this work was to demonstrate further the 
applicability of the FEM to the non-linear reactor dynamics 
problem a3 well as toinvestigate the dynamic response of the 


reactor fuel subassembly. The analysis required a novel 
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approach to handle the gap conductances present at the in- 


terfsces of the equivalent cell model of the fue] Subasseme- 


bly; this will be clarified in the analysis. 


Ti, OESCRIPTION OF PROBLEM 


A. PHYSICAL SYSTEH 
The typical Liquid Metal Fast Breeder Reactor (LMFBR) 


‘' 
4 


core consists of many hexagonal modules, each containing 


Seyeral hundred fuel pins. For this analysis, an equivalent 
cylindrical ceil is used to mode? the fuai subassembly: see 
Figure 1. The use of equivalent cells as models for larger 
Systems has been common practice in nuclear analysis (4.e., 
the well known: Wigner-Sietz method). In using an equivalent 
cell, the actual shape of the reactor core is not important, 
and the analysis is applicable to any reactor which has the 
= same equivalent cell. 

The equivalent cell considered in this analysis, Figure 1, 
is fueled with enriched uranium dioxide, has a stainless 
steel cladding, and has liquid sodium for a coolant. The 
dimensions used are 

a= 0.254 cm, 

b = 0,292 cm, 

c = 0.365 cm, 
and H = 33.0 cam, 
The gap between the fuel and cladding is very small and, in 
fact, may be nonexistent as in bonded fuels. The dimension 
of this gap has been assumed negligible, The height, H, of 


the fuel rod is shorter than many proposed systems (Fast Flux 


Testing Facility and Clinch River Breeder Reactor). However, 


COOLANT FLOW 


CLAD-COOLANT 
J INTERFACE 
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Figure 1, Equivalent Cyliadrical Cel] 
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ta facilitate the numerical solution a smaller rod was used, 
The dynamic behavior prior to fue: 2's failure for this sys- 
tem should be similar to the Ae se <v =F larger systems. 

The treatment of this problem in three dimensions would 
be prohibitive in computer usage. “Zerefore, azimuthal sym- 
metry is assumed, and the problex cscomes a two-dimensional 


cylindrical (r,z) problem. 


8. SYSTEM MODEL 
The analysis considers the monoénoargetic neutron diffu- 


ston approximation to model the transient neutron transport 


problem. A simple conduction-convection heat traiis¢ar model 


is used for the energy transport probiem., The temperatures 
in the model are direct*; coupled to the neutron flux througk 
the nuclear heat generation within the fuel. The neutron 
population {s in turn coupled to the tamperature through any 
of a number of reactivity feedback mechanisms. The nuclear 
Doppler effect is perhaps the most important of these mecha- 
nisms since it provides a negative temperature coefficient 
which increases the inherent stability of the reactor. Prior 
to core disassembly and fuel melting, the nuclear Doppler 
effect is the most dominant feedback and is, therefore, the 
only feedback mechanism consideved in this analysis. 

For irradiated, mixed-oxide fuels, a phenomenon of fuel 
restructuring has bees commonly osserved. This restructuring, 
essentially a@ change of phase of the fuel nzterial, presents 
a unique heat transfer problem particularly during transient 


conditions. The problem has not been fully characterized and 
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is beyond the score of this analysis. The fuel nae: there- 
fore, assumed to be a homogenevus mixture cf er:ichked uranium 
dioxide, 

At the fuel~ciadding interface, there exists a gsp which 
produces a therméi resistance. This thermal resistanca is 
one o? the most signifscant deterrents to the unergy transfar 
to the coolant. The interface may be in physizal contact cr 
an actual gap may exist. The prediction of the thermal resis- 
tance is extremely complicated and must take into considera= 
tion many parameters: initial dimensions, type of bond, fi11 
gas compositior, fuel restructuring, fuel swelling, prior fuel 
life cycle, toa mention a few. Reference [5] documents a com- 
puter program which attempts to precict the gap canductance, 
Noap* 
the same manner as a convection hsat transfer coefficient when 


This treats the thermal resistance at the interface in 


considering convection heat transfer, Since it is not the 
objeztive of this analysis to predict the gap coefficient, a 
representative set of vaiues for gap coefficient, as given in 
Ref [6j, is used in this analysis, These vajues ars assumed 
to remain static during the transient. The gap conductance 
profile actuatly varies with time and will have an effect uoon 
the transient, &3 noted in Ref. [7], The prediction of this 
variance was no‘, considered important for this analysis; 
therefore, the static assumption was made, 

An average convection heat transfer coefficient was used 
to determine the heat transfer from the ciadding to the cool-~ 


ant. The value used was determined from an empirical formula 


given in Ref. [5] and repeated in Annendix C. 


In this work, consideration is given to step and ramp in- 


creases in reactivity, aithough any reactivity transient may 
aasily be considered. The step and ramp increases in reactiv- 


ity probably represent the most reatistic physical reactivity 
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inputs in a reactor. Once the reactivity has been inserted, 
the transient overpower excursion begins. Unless the Doppler 
feedback can override the inserted reactivity, the excursion 


will continue until there is physical core disassembly. 


C. NUMERICAL SOLUTION 
The system of equaticns which models the proposed problem 
is a non-linear, initial-boundary-value problem. The conven- 
tional method of solution of the reacto: dynamics is the point 
kinetics formulation. It was pointed out in Refs. [1], [2], 
[3], and [4], that there is a non-negligible error in this 
model, particularly under conditions of asymmetric disturbances 
or space-dependent feedback, Reference [4] demonstrates the 
somewhat novel approach of using the finite element method 
(FEM) to solve the space-time dapendent reactor dynamics prob- 
jem. As shown in Ref. [4], the FEM is quite effective in 
handling localized perturbations and space-dependent feedback. 
in this work, only uniform disturbances were considered; how- 
eves. the feedback model was space-dependent. Therefore, the 
finite element method is used to solve the non-linear, coupled, 
space-time dependent neutronic and heat transport field equa- 
tions. The solution technique results in a large computer 


storage requirement; therefore, an optimum compact storage 


Scheme, Ref, [8], is utilized for storage of the discretized 


matrices. 
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Once the domain has been discretized by the FEM, the solu- 
tion of the resulting ordinary differential equations was to 


= be accomplished by the implicit Gear's method, Ref. [9]. 
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Tif. MODEL DEVELOPMENT 


A. WEUTRONIC ANALYSIS. 

In this section the governing field equations for the 
neutron population (flux) for each of the three regions 
(fuel, clad, and coolant) of the domain will be formulated. 
The monoenergetic, diffusion theory will be used. 

Consider an arbitrary volume of material within a reac- 
tor. Applying the condition of conservation to the mono- 
energetic neutrons leads to the neutron equation of 
continuity [10] 


an(r.t) = S(r,t) - B.(r)e(rst) - div d(r,t) (1) 


where 


r= spatial point 
t - time 
n(r,t) - neutron density 
S(r,t) - neutron production 
z,(f) - neutron absorption cross section / 
o(r,t) - neutron flux 


J(r,t) - neutron current 


The left-hand side of equation (1) represents the time rate 
of change of the neutron density which is related to flux, 
9 , by 


n(rt) = + 9(r,t) (2) 
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where 


Vv - neutron velocity 


On the right-hand side of equation (1), the first term is 


neutron production, the second term is a neutron loss through 


absorption, and the third term is a neutron loss through leak- 


age from the control volume. 


Using equéti.1 (2) and applying Fick's Law to the equation 


of continuity results in the classical neutron diffusion equa- 


tion 
ve(D(r)vs(rst)) -E,(rio(rst)+s(r,t) = 5 Sgltot) (3) 
where D(r) - neutron diffusion coefficient 


The vector notation used here is intended to include only two 


dimensions (r,z) since azimuthal symmetry has been assumed. 


Equation (3) {s applicable to each of the three regions 


of the equivalent cylindrical cell. Zhe subscripts F (fuel}, 


c (cladding), and co (coolant), will be used to denote these 


regions, 
1, Fuet Region 
In applying equation (3) to the fuel region the 
material properties of the fuel must be used. Within the 
fuel the neutron souce term is due to the nuclear fission 
Process, DOuring fission, neutrens are released as both 


prompt neutrons and delayed neutrons so that 


S(r,t) * So (ret) o Sp(r.t) (4) 
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where 55 (n>t) ~ prompt neutron source 
Sp{r.t) ~ delayed neutron source iq 


The neutron sources are commonly represented as [10] 


Sy" he Ege Op(1+8) (5) 
and 
n 
Sy "4 Cia, (6) 
where 
k, - infinite multiplication factor 
B - fraction of fission neutrons which appeer 
as delayed neutrons 
nh = number of delayed neutron groups 
C; ~ Concentration of delayed neutron precursors 
in the {th grous 
A, ~ decay constant of the delayed neutrons 
The space sad time variables, r and t, wil) be Cropped except 7 


Where needed for clarification. 


The concentration of delayed neutron precursors, Cys 


is given by the following first order partial differentia] 
equation [10] 


aG yi} 
TE = Bake Eapop ~ Ayly (7) 


where 


8; ~ fraction of delayed neutrons which appear 7 
as defayed neutrons in the ith group 


The solution of equation (7) is 


ee ; 
Cy By f MEEK Crea pepctte cre t (8) 
0 


ce - concentration of delayed neutron precursors 
of the ith group at time zero 


Reference [10] develops an expression for the initial concen- 


tration of delayed neutrons 
CP = B KY Zap O°/A, (9) 
where 4 


kS - initial finite multiplication factor 
6° - initial neutron flux 


Combining equations (3), (4), (5), (6), (8), and (9), yields 


the governing equation for the fuel region 


Vo (O-%¢e) + Eapoplk(1-8) <1] 


n 
+ Dajte, fe Ag(EtN), (et ogd, pdt! 
i=] 0 
op 
“i,t 1 °TF 
+ Bike Eapepfdy OUT * 7 aE Sa 
2. Cladding Region 
The cladding separates the fuel and coolant and con- 
tains the fuel and the fission by-products. Equation (3) 


governs the neutron flux in the clad, Since the clad contains 


no fissile material, the neutron source term is zero. The 


field equation is, then, 


< 1 3% 
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3, Coolant Region 
The annular region around the cladding of the equiva- 
Tent cell contains the coolant. As in the clad, there is no 
fissile material in the coolant and, therefore, no neutron 
source term, The equation governing the neutron flux in the 


coolant is 


, 96 
V*(DegVbeq) ~ Lacotco * ; at - M2) 


Equations (10), (11), and (12) are the one-velocity, 
diffusion approximation used to model the neutron transport 
problem. 

4, Infinite Multiplication Factor 

The infinite multiplication factor, k, may be ex-~ 
pressed as the infinite multiplication factor at time zero 
(start of transient), ko, plus the postulated reactivity in- 
sertion (such as a step or a ramp), p, minus the change in 
the Doppler reactivity feedback, App. Other feedback mechan~ 
isms are normally not as significant as the Doppler broaden- 
ing feedback prior to fuel melting and have been neglected 


in this analysis. Therefore, 


k, * ky + 9 = App 


For a fast reactor, kS may be approximated as 


mM 
> 
n 


ke xy 


“4 
i] 
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v - average number of neutrons released 
per fission 


Leg - fission cross section of the fuel 


rae - absorption cross section of the fuel 


The nuclear Doppler effect is a very important safety 
feature in @ nuclear reactor, Nuclei in an atom are in con- 


tinual motion due to their own thermal energy, As a result 


- 


of this motion, even when monoenergetic neutrons interact 
with the atom, there appears to be a spread in the energy of 


the neutron ~ the Doppler effect. It can be shown that the 
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cross section of a resonance becomes less in magnitude and 


wider as the motion of the nuclei increases [10]. As the 
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temperature increases, the motion increases, and the shape 

of a resonance cross section broadens, This broadening in-~ 
craases the average cross section, thus, providing a negative 
temperature coefficient. This effect is shown in Figure 2. 
It is this nuclear Doppler broadening effect which provides 
one of the few inherent, reliable, negative reactivity feed- 


backs which slows an overpower transtent and possibly stops 


a mild overpower transient. 
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Figure 2, Doppler Broadening of a Resonance Peak 
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The Doppler reactivity change with respect to fuel 


temperature changes should be written as [6]. 


ss wat 3/2 ¢ ptm) + ttl (15) 
where a, b, and c are parameters determined from experi- 
mental work and m is an integer. However, as noted in 
Ref. "6], a substantial amount of work has shown that 

T a is very nearly constant over the temperature range 
under consideration. Therefore, vne coefficients a and 
c have been set equal to zero, and b is defined as 


dp 
b= ky * Tae (16) 


The constant, Kos is commonly called the Doppler constant. 
Solving equation (16) for Py yields 


where K is an integration which may be obtained. from 
initial conditions. 


Kk = pg - b In Te (18) 


where 


5 - Doppler effect at time zero 


Te - fuel temperature at time zero - 
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Substituting for K in equation (17) will give 


Py - Py = Apy = b In(Te/Te) (19) 


The infinite multiplication factor now becomes 


k, = k& +p = b In(Tp/T2) (20) 


The effect of delayed neutrons is small compared to 
the prompt neutron effect; therefore, it may be assumed that 
the Doppler effect on delayed neutrons is insignificant to 
the overall problem. The k,. of equation (8) is, then, 
assumed to be k° , With this assumption and equation (20), 
the neutron diffusion equation in the fuel, equation (10), 


may be rewritten as 


Ve(DeVor) + Eade lka(t-8)-1+p(1-8)-b(1-8)1n(Te/TE) ] 
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(21) 


To facilitate the prasent analysis, the number of | 
dalayed neutron groups is taken as one averaged group. So 
that, A; and 6, become \ and B, respectively, This approxi- 
mation should have little effect on the problem under con- 
Sideration. 
5. Boundary Conditions 
The neutron diffusion problem involves solution of 


the partial differential equations (11), (12), and £21), with 
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the following boundary, interface, and initial cenditions: 
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1) se (0,z,t) = 0 i 
2) op (a,2,t) = o, (asz,t) 
3) Dp bd (a,z,t) = D, ne (a,z,t) 


4) %e (b,z,t) = dog (bsz,t) 


3h ag 
5) 0, ape (b,z,¢) = 0, x S2 (bz, t) 


r] 
6) g8%e,2,t) = 0 
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7) plot Btls octrat Ht) = dcg(rat Het) = 0 
| 8) bp (10) = 42 (x) 
9) , (n.0) = 92 (r) 
10) oo(ns0) = $25 (r) 


Boundary conditisn 1) results from the assumed azimuthal sym- 
metry. Interface conditions 2), 3), 4), and 5) are continu- 
ity conditions of the flux. Boundary condition 6) results 
from the use of an equivalent cel] and basically indicates 


there is an equal number of neutrons transferred in and out 
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of the cell at the outer boundary, This should be valid un- 
Tess the cell is located near the outer edge of the reactor. 
Boundary condition 7) is an assumption that the flux is zero 


at the axial boundaries of the cell. Initial conditigns 8), 


9), and 30) are the assumed initial distributions of the 
neutron flux, 


B. HEAT TRANSFER ANALYSIS 


In this Section, the principle of conservation of energy 
will be used to formulate the governing ftald equations for 
the heat transport in each of the three regions. A simple 
heat conduction model with convection heat transfer to the 
coolant is used to model the heat transport problem, A gap 
conductance mode] is used to describe the heat transport 
across the gap at the fuel-clad interface, 

1, Fuel Region 

Conservation of energy within the fuel region yields 


the unsteady heat conduction equation with a generation term 


oT 
Ve(ke(r) OTe (ret) )+a(ret)*oe(r)Cye(r) gee (ret) — (22) 
where 


ke (nr) - thermal conductivity of the fuel 
Te(n,t) - fuel temperature 


4(r,t) - nuclear energy generation per unit 
volume . 


Pe(r) fuel density 


Cop (n) - fuel specific heat 


As in tha neucronic analysis, the vector notation is intended 
to include only two dimensions, (r,z). The r and t will be 
dropped except where needed for ciarificaticn. 


Tne nuclear generation term may be expressed as 


Ge Lee of (23) 
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where @ - nuclear energy released per fission 


As can be seen, it is through the nuclear generation term 
that the temperature is directly coupled to the neutron 
flux. This. coupling and the temperature dependent Doppier 
reactivity feedback combine to make the coupled problem non- 
linear, 

Substituting equation (23) into equation (22) yields 
the governing thermal equation for the fuel 


aT 
: F 
Ve(keVTe) + @ Eee Op * OF Coe Em (24) 


2. Cladding Region 

Conservation of energy within the clad will yield 
the heat conduction equction. With nuclear generation, a 
relatively small amount (~5%) of the energy will ce released 
in the cladding and the coolant. However, in this analysis 
it is assumed thet the total energy release is in the fus! 
region, This should rot create any significant error. Us- 
ing this assumption, the unsteady heat conduction equation 
for the cladding becomes 


| aT, ; 
Ve(kKVT.) = be Coc BE (25) 


3% Coolant Region 
Once again, conservation of energy will lead to 


the heat conduction equation plus an additional term which 
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takes into consideration the coolant flow. 
equation is 


The governing 
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V*(keg¥Toq) - Veo Fz! coeae0% co) co “nco FE 
(26) 


where Veg 7 Coolant flow velocity 
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Equations (24), (25), and (26) are the governing 
equations used to model the energy transport problem. 
§, Interface Conditions 
The interface between the fuel and the cladding may 
be an actual gap with a finite distanc:, or the surfaces 


may be in intermittent contact on a microscopic scale. To 
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model the hest transfer e¢coss this interface, a tp heat 


transfer coefficient ¢2 introduced. The gap cosfficient 


must take into consideration many items (e.9., radiation 
heat transfer across the gap, heat transfer %y solid-to-solid 


contact, heat conduction across a gas filled gap). The pre- 


diction of this gap coefficient is extremely complicated and 


beyond the scope of this work. In Ref. (6], @ set of values 


p is given and the axial variation of Heap in this 
analysis is approximately the same. 


for Hoa 


A cosine curve has been 
fitted to the sample data to determine the gap coefficient, 


see Figure 3. The heat flux across the fuel-clad interface 


is, then, 
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Q* Hay (Tp - T) (27) 
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The heat conducted out of the fuel is governed by Fourier's 


equation and is equa? to the heat transferred across the gap 
aT = 
q = oke a 28) 
Equating equations (27) and (28) gives the fuel-clad inter- 
face condition 


ke(a a,z) ate 
Te(a,z,t) + He, p te) ala, z,t) = T (a,z,%t) (29) 


and from continuity 
ate oA. ee 
ke(a,z) ap-(az,t) . k(asz) =x7(4,Z,t) (29a) 


As is common practice in convection heat transfer 
analysis, 2 coolant surface heat transfer coefficient, hsurf, 
is used to account for the thermal resistance at the clad- 
coolant interface, Similar to the fuel-clad interface analy- 
sis, the clad-coolant interface conditions may be determined 


k.(b»2) wr 
Telbozst) + FE grElbyzyt) «Too(bszet) (30) 


aT. aT KA 
k.(b,z) ap7(b»z,t) = keg (bsz) Br (bez t) (30a) 


5, Boundary Conditions 
The boundary and initial conditions for the heat 


transport problem are: 


sm ult 


ri 
2) Tyg lte- got) = Tor enuy 


T 
3) it (r, x,t) = 0 


Fy 
bell 
ae 
“| <4 
Oo 
oS 
— 
is] 
ww 
nN 
ww 
ce 
Lael 
a 
ren) 


5) = (r,+ Bt) » Te (r+ St) = 0 
gz od oe gz (ME Bs 
Te(r) 

Te(r) 


TQ (n) 


ba I fo) 
al baal 
— =< 
ir] ba 2 | 

tt ~~ 
[x Iss 
= . 

a o 
baal ~~ 
” " 


ie 4 
ee 
~d 
~ 
be 3 
~~» 

=] 
atl 


Boundary condition 1) results from the assumed azimuthal sym- 
metry. Tha coolant has been assumed to enter the flow 
channel at a constant temperature, TpLENUM’ This results in 
condition 2). Soundary conditions 3) and 5) result from an 
assumption that no heat is transferred axially from the fui's 
rod, Boundary condition 4) is the result of the use of the 
equivalent cell. Conditions 6), 7), and &) are the assumed 
initial conditions. 

Solution of the nonlinear coupled neutronic and 
energy transport problem involves the solution cf the partial 
differentia] equations (11), (12), (21), (24), (25), and 
(26), with the appropriate boundary, interface, and initial 
conditions. 

Several works, Refs. [4], [11], [12], have demon- 
Strated the feasibility and success of the finite element 


Method in solving nuclear reactor dynamics problems, The FEM 


is used to reduce the partial differentiai equations developed 
in this analysis to a system of ordinary differential equa- 


tions. Integration of these ordinary differential equations 


{ODE} yields the solution. 
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IV. FINITE ELEMENT FORMULATION 


In this section, the basic theory underlying tne finite 
element method is formulated. Selection of the finite ele- 
‘ments and the shape functions for the elements. are given. 
Some simple transformations which facilitate the integration 


necessary in the FEM are also presented, 


A. BASIC THEORY 
To obtain a numerical solution, the governing partial dif- 


ferential field equations are transformed into a system of 


ODE in finite dimensional vector space. This may be accom- 


plished in several manners such as the finite-difference 
method, the variational method, or the weighted residual 
method. In this work, the Galerkin method (a weighted re- 


sidua] method) is utilized for the discretization of the 
spatial domain. The Galerkin procedure will be applied to 
each of the governing field equations. Any of these equa- 


tions may be considered to be in the following form 
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Sy (rst) - Ly(rst) = f (ret) (31) 


where ww represents the unknown function, e.g., in equation 


(21), ww represents OF » £ represents the operator for 


each individual equation, and f is a forcing function. In 


the finite element method the solution is approximated as 


N 
y(r,t) a h(r,t) wes Ny (1) y(t) = <Ny>thy} (32) 


where N = the number of degrees of freedom 
Ny - the »lament shape function 
vy unknown coordinate function 
< > ~ matrix notation for a row vector 
{ } - matrix notation for a column vector 


<Ny> bs <N, No oes Ny eee Ny> 


The residual, R(r,t), is a measure of the error in this 


finite element approximation. The residual may be considered. 
as 


R=ae-fv-f (33) 


The best solution for % is one which "minimizes" this 
residual. Various "minimums" are obtained by the weighted 
residual methed by setting 


fale) RAV =O HT,2y 000K (34) 
V W,(r) - weighting functions 
With the Galerkin method, the weighting functions are the 
shape functions defining the approximation of equation (32) 
(i,e., W, = Ny). A noteworthy attribute of the Galerkin 
method is the opportunity of using an integration-by-parts 


anasto a ite tl i sl) i it Us ota saps 
sant asa nes iS ns a gt all ela : 


of the terms involving the second order spatial derivatives. 


A Tower order finite element may be used than would have 
been possible otherwise. 


been chosen, the problem becomes 


J N, (Se -fy-f} v= 0  i#1,2,...8 (35) 


The integration involved in equation (35) is carried out 


on the element level, taking advantage of the use of a “local” 


coordinate system. Once the integration is accomplished, 


the results are merged into a system using "global" coordi- 
nates. On the element level 


z° = <Nj>® 953" JET pe pean (36) 


where the superscript e incicates the element level. Sub- 
stituting ¥° into equation (35) and noting Cy 53° ig not 
a function of the spatial domain yields 


Sera [<ny>® £(Ns}®~<Ns>*#*]ave(psd® = 0 
(37) 
where ae | -” Vi2sevall® 


N® = number of degrees of freedom for 
an element 


The operator J will vary depending upon which governing 
equation is under consideration. 


Once the weighting functions have 
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: B. SHAPE FUNCTIONS 
The shape functions, Ny» are chosen to satisfy certain 
A completeness and convergence criteria [13] and will depend 
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upon the finite-element used for the spatial discretization. 
iF 

Many previous works, Refs. [4], [11], and [12], utilized u 

| linear triangular shaped elements to discretize the spatial 


domain, This. element was the first element considered. | 


However, because the width of the cladding is very thin, 
elements in the cladding region would have extremely large 
aspect ratios (ratio of base to height) unless an petrenely 
large number of elements in the axial direction were: used, 
A large number of elements becomes numerically untractable. 
7 Previous experience with triangular elements had shown that 
| large aspect ratios yield inaccurate resuits. Ziamal, Ref. 
[14], showed the error, e, when using triangular elements, 
1s proportional to the square of the longest side, h, and 


inversely proportional to the sine of the smallest angle, y 
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A triangular element with a large aspect ratio necessarily 
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must have a small related angle which adversely affects the 
error in the FEM. Hopefully to alleviate the problem, an 

fsoparametric, quadratic, rectangular element was selected, 
The aspect ratio would still be large, but experience, Ref. 
[15], with the use of rectangular elements indicated that a 


large aspect ratio is not always a detrimental factor, 


The shape fu ctions for this element are wel] documented, 
Ref. [13]. utilizing a “lceal" coordinate system (See Figure 
4), the shape functions may be written as 


Corner nodes 4 = 1,3,5,7 


Ny * #(14E,) (14n,) (Etng@2) (38) 


Midside nodes Ny * (1-62) (14n,) 1=2,6 (38a) 
Ny = (14E,) (Ton?) 194,8 (38) 


where Bo = gS, 


These normalized shape functions are shown in Figure 5. 


The jocal and global coordinates are related by the 
following. 
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r= <Ny>® (r,3® (39) 
and z* <N,>° {z,}° (39a) 


C. COORDINATE TRANSFORMATIONS 
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When using a Jocal coordinate system, some simple trans- 
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formations facilitate the integrations required by equation 
(37). In cylindrica? coordinates with azimuthal symmetry 
dV = 2nrdrdz (40) 
The derivative terms may be transformed by the following d 
- cay") (41) { 
any j 
Le F 


GLOBAL ELEMENT 


Ue g 


LOCAL ELEMENT 


AFTER TRANSFO/MATION 


ete, 


Ya~] 
LOCAL NUMBERING 
SYSTEM 


Figure 4, Elament Transformation 
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Figure 5. 
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where fa}! js the inverse of the 2x2 dacobian matrix de- 
fined in Appendix A, A$ shown in Appendix A, this inverse 
can be easily shown to be (for this problem) 


* * 

Jy, 412 

tay) (42) 
* * 
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where 
* : * 
U1 s 2/ar r Yi2 x 0 
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Jo} = 0 : Joo ba 2/Aé 
radial length of the element 
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AZ = axial length of the element 
Elements of area transform as 
drdz = det[J] d&dn {43) 


For this particular problem, det{J] may be shown to be 
(Appendix A) 


dettay = A (44) 


where A® ~ area of the element 


Elements of sxial length become 


; Le 
dz = 5- (45) 
; where L® ~ axial length of the element 
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Utilization of these transformations makes the integrations 
required fn equation (37) amenable 


to integration by numeri- 
cal Gaussian quadrature, 
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¥. APPLICATION OF FEM TO GOVERNING 
FIELO EQUATIONS 


In this chapter, equation (37) is applied to the govern- 
ing field equations previously derived. Tha eiement matrices 
for each of the operators are developed so that the discreti- 
zation of the spatial domain may be accomplished. The inte- 
grations required by the application of equation (37) are 


performed numerically using Gaussian quadrature. 


A. GAUSSIAN QUADRATURE 


Prior to the application of equation (37), it is appropri- 
ate to discuss briefly the procedure used for the numerical 


integration. The product. Gaussian quadrature formsts is [16] 


iP. ; m OM 
I, * Sf 9(E,n)dndé = 5° > WyWs 9(&4 ons) -(46) 
-1 -1 f=] jel 
where I, - area integration 
g(E,n) - any function of — and n 
Wy 4. -* weight associated with location 
: for j 
a - number of Gauss sampling points 


in one-dimension 


The values of the weights associated with aach Gauss point 
are given in Ref. [16]. Equation (42) may be simplified 

somewhat by combining the summations and weights 
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For line integrations, the Gaussian quadrature formula in- 


volves only one summation 


a 1 s 
1 J f(n)dn 7 ie (nz) 


B. NEUTRONIC FIELD EQUATIONS 
The discretization of the spatial domain by the finite 


element method is accomplished by applying equation (37) to 
the governing field equations, using 


oe * <n >* (¥,5)° J21 42500038 
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1. Fuel Region 
The govern’iig equation for the fuel region, equation 


(20) after applying squation (37) becomes 


N, 3¢ a6 
anf fit = rdrdz-2x f f Ni {% Fp (rp art) 
rz rz 
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(50) 
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The second order terms in equation (50) may be re- 


duced to a first order term by application of Green's Theorem 
or equivalently integration-by-parts (See Appendix B). 
Dividing through by 2m and reducing the second order terms 
yields 


26, 7 Ad 3 ma Wd 2 
f rNyDe aoe dz + frnyde aot dr [fz spr rdrdz 
r rz 


aN, 2p ONY 36 
hf {Petar ore tage aged Eg php( 1-8) [kgto-bin(Tp/Fey)] 


t 
AyKBEyy JME) (Rata dygae ayer ht rdrdz = 0 (51) 
0 


From continuity and boundary conditions the line in- 
tegrals are zero, Now using the approximate functions of 
equation (49) and noting that {yp)® is not a function of 


space, equation (51) may be written as 
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Se reugreratep" - Xere7t ff {N, }rdrdz = 9 (52) 
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The integrations may not ba easily carried out with 


a shift to local coordinates and with use of the previously 


derived transformations. Rearranging equation (52) and 


assuming the properties are constant for each time step gives 


ia & x * * d e 
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= td 


if ip et pprdetLd]eedn p}* = 0 (53) 


Since the element chosen has eight degrees of freedom 


(nodal points), the discretized matrices which result from 
the integration of equation (53) will be 8x8 matrices and the 
forcing function will be an 8x1 vector at the element level, 
Defining the matrices as 

J * * * + i 
a a CONy Jy) }<34, Ny er tthy Jae }<Joo Ny n> srdetlo] d&dn 
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HME DCN nde 922 Fey (54) 


4 
a 
a 
q 
a 
5 
i 
3 
j 
a 
a 
a 
; 
4 


§ a 8 
where u, » HW and ry pie ri(Hy), 
(= 


17 ¥ me 


2 
11 ae as 
ae {Hy Irdet(J]dédn = {Fi}g.q * 7 P CALE (56) 
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J f In(Tp/Tp) {Hy HN ,ordet[I]dgdn * [H4, 5], 
Ate. i 4 : 
= § Loire Teed HM Tae (57) 
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To carry out the summation of equation (57}, the 
temperature Te must be known; however, th? temperature is 
exactly what is being sought, To alleviate this probles, 2 
linearization is used. 

In the solution technique, the temperature {is pre- 
dicted for the next time step. It is this temperature which 
is used for the determination of matrix H4,. 

Equation (53) simplifies to 


- 


(DpLH12, ,T+L2ye(1~8) (kote) XBL, f(t) ICH, #2, -(T-ADBIHG, sTHGp* 


“lcve"*t ry + Mua, 1G} = 0 (58) 


The function f(t) is evaluated by summine the 
values at each time step using the trapezoid rule for numeri- 
cal integration. 


t «© 
f(t) = evAt y ett Tet ddt! 
0 


Defining 
1,Co(t)} = 5 hyto(ty)+ 9(%4_4)) 
g(t) = eb t(Ke+p) 
hy time step taken 
The function may be expressed as 
s 
f(t) map» 1, [g(t)] 


S = number of time steps 
2. Clad Region 
Following the same procedure as with the fuel equa- 


tion, the discretized form of equation (11) becomes 
{D,(H12)5] + EaclH3, sTHy_}® + FHS, sIGH,}® = 0 (62) 


3. Coolant Region 
The governing equation for the coolant region, equa- 


tion (12), may be discretized into the form 


{D,gfH12, 4] i" Bacolt3s j]}H¥eqh"+ HHS, j1ligg)”* 0 (63) 
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Once the governing equations have been discretized 
at the element level, they are combined into a system of 
equations, at the global level. On the global level the 


governing equation for neutron transport takes the form of 
CHnxnf¥nxy * CPInxnf¥inxt * fFiayy * (64) 


where [HJ], [PJ], and (FJ represent the system matrices and 
n is the number of nodal points used in the discretization, 
There are, then, n simultaneous ordinary differential equa- 


tions used to describe the neutron transport problem. 


C, HEAT TRANSPORT FIELD EQUATIONS 


The spatial domain for the heat transport problem is 


discretized in the same manner as the domain for the neutronics 


¥Yelid. Let 


T . a J21,2,.+.58 (65) 


problem. The same element matrices pre'iously defined are 
k=F,c,co 


1. Fuel Region 
The governing field equation for the fuel, equation 
(24), is discretized by applying equation (37). Using an 
integration by part to lower the order of the second order 
terms allows equation (24) to be written as 
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From continuity considerations the line integrals 


are zero #xcept along the boundaries of & region. It is 
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assumed that no heat is transferred from the cell in the - 


( 
2 
~ 


axial direction (boundary condition 5); therefore, the 
first line integral of equation (66) is zero. In the neu- 
tron diffusion problem there was continuity of flux at the 
interfaces so that the line integrals were zero; however, 
the heat transfer at the interfaces is affected by the gap 


and film conductances, The fuel-clad interface condition, 


ate 
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equation (29), may be rewritten as 
oT, aT. 
oke gp * Hgap(TeeTe) = ke a (67) 


Substituting equation (67) into (66), dividing by minus one 
and utilizing equation (65) yields 
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Transforming to local coordinates and integrating 
by Gaussian quadrature yields the same element. matrices as 
given for the neutron flux, equations (54) and (55), except 
for the line intesrals between regions. It should be noted 
that the line integrals exist only on the interfaces, along 
which there is a disconZinuity of temperature. For the 
fuel equation the interface corresponds to the local co- 
ordinate &=1. Define 

2 


1 en 
{ Wyletipdn = Okt; seg * Z »» (Ny), (Hy) My (69) 


where the N's are evaluated at &=1, Many of the terms of 

Kl will be zero since only the nodes on the €#1 boundary 
will have shape functions which are non-zero, It is through 
Kl that the temperatures for each region are coupled to- 
gether. Equation (68) may now be written as 


Pat gap CERI (tp}®-CK1 Jr }O D+keCHI2] (6) 
weLecLH3](ve}® + opC LAB] {te) = 0 (70) 


In obtaining this equation, it was assumed that material pro- 
perties for each nodal point were constant at each time step, 
Perhaps a better assumption would have been to assume an 
average value for the properties of each element. The dif- 


ference should not be significant, and the assumed constant 


nodal properties were numerically more tractable, 
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2. Clad Region 
Apslying equation (37) to the governing field equa- 


tion for the clad, equation (25) gives the discretized form 
of the equation. Assuming no haat. transfer in the axial 
direction on the boundaries (Boundary condition 5), equation 
(25) becomes 


Te aN, aT, aN, aT, 
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aT. 
Fe Soe Ny ape tdrdz = 0 (71) 


In the cladding, there are two interfaces along which 
the line integral of equation (71) {is not zero, along the 
fuel-clad interface and along tne clad-coolant interface, 

For the fuel-clad interface, equation (61) applies. For the 
Clad-coolant interface, the interface condition, equation 
,30), may be rewritten as 


aT. aT Ko 
Ke Tr =» hsurf (ToT) = wk + (72) 


Along the fuel-clad interface, the local coordinate 


corresponds to —=-]. Define the new element matrix 
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where the N's are evaluated at &=-1. As with K1, K2 will 


have fiany Zero values because the shape functions are 


evaluated st é=-1, 


Aiong the clad-coolant interface, the local coordi- 
nate corresponds to €=1 and the Ki matrix is appropriate, 

Substituting equations (67) and (72), equation (71) 
becomes. 
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The governing equation for the clad region may now be written 
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pe 
The line integrals of equation (73) affect only nodes which 
are on one of the boundaries; therefore, the nodal inouts 
into Kl and K2 are zero unless the node is on one of the 
boundaries, 
3. Coolant Region 
The field equation governing the coolant may be 

discretized in the same manner as above. After applying 


the Galerkin method and performing an integration by parts 


on the second order terms, equation (25) 
becomes 


C on, aT, aH, aT 
a co rl, dz + i 2 Kool at rue! rdrdz 
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The lina integral, when evaluated atc, is zero (boundary 


condition 4). When evaluated at b, or correspondingly at 


E=-1, equation (72) is valid and K2 matrix is appropriate. 


All the terms of equation (75) have been defined except the 
flow term. Define 
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Transforming to local coordinates and integrating raduces 
equation (75) to 


MpMeu ng (KZ tog) [KZ tH") + kK, LHIZ]t,,}° 


é a ud e. 
* VeoPcompcottSlitg}” + Peopeolts]{t gh = 0 


rl au ds Wc a ie li Sica Ai ge sa a a sho lp gba taiaibe, ti 
1 ost glut "Latitogl ot Sugg halal dl oil Pee Tee ee i ls ra 
sige thnks YY MAD suis 
if a a i tl at Sul Shao 


dale hi ALS Wy tga 


Now that the governing equations have been defined 
for each region on the element level, equations (70), (74), 
and (77), they may be assembled into a system equation on 
2 the global level. The equation will be in the general form 
of 


CK ayn lthaxt . CMI xn ¥ axl . CG Yaga fthaxy a 


(78) 


D. DISCRETIZATION CF THE SPATIAL DOMAIN 
Prior to the numerical solution of the governing equa- 
tions, equations (64) and (78), the spatial domain must be 
divided into a number of elements. ‘or this work the domain 
was subdivided as shown in figure 5. 
Since there is a discontinuity of temperatures at the 
interfaces, as described by equations (29) and (30), a 
noval application of the FEM method was necessary. The com- 
mon practice for handling these "flux" type boundary condi- 
tions is tu define a cunstant reference temperature, T, ; 
as whan working with a convection heat transfer problem [17], 


or to define a known function, as when working with a 


fracture mechanics problem [18]. In either case the refer- 
ence condition was known. The novel appiication here lies 
in the use of a different field equation to describe the 
reference temperature, e.g., the clad equation (74) is the 
reference condition for heat transfer from the fuel across 


the gap interface. 
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The discontinuity of temperature at the interface neces- 
sitated another novel application of the FEM. Since there 
is a temperature drop aleng each interface, 2 single node 
there is not adequate. In the discretization of the domain, 
two nodes were used for each interface point (for example, 
points 62 and 63 in figure 6). This aliows the temperature 
drop due to the gap and film conductances to be taken into 
consideration. Since two nodes are used, the governing equa- 
tions for each region are not directly coupled together. 

The coupling of the regions is accomplished by the "flux" 
boundary or interface conditions since it is assumed that 
any heat flux leaving a region enters the adjacent region. 


Consider a typical set of elements on an interface 
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LAB-COOLANT INTERFACE 
FUEL-CLAD INTERFACE 


Figure 6, Finite Element Discretization 


The coupling terms Ki and K2 may be combined into a system 
K matrix which shows the coupiing. The K matrix for the 


Simple set shown is 
12345 6 7 Bee 13 4 15 6 
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As can be seen, the nodes con the interface are coupled to 


the adjacent element interface nodes. For example, node Z 


in element I is coupled to nodes 3, &, and %4 in element II. 


E. OPTIMUM COMPACTING SCHEME 

tc.) are nxn matricas, 
where n {#3 the number of nodal points used in the discre- 
tization of the domain. In terms of computer storage, these 
matrices may become excessively large if they are stored as 
nxn. There are several techniques avajlable to reduce this 
storage requirement. The most common method {is the banded 
storage scheme, whereby only the banded portion of the 


matrices are stored. With judicious numbering of the nodes, 


considerable savings may ba realized. Howaver, it is not 
the optimum storage scheme [8]. 

Since the shape functions, N,, for the pth nodal equa- 
tion are nonzero over only the element containing k, the 
System matrices are not only banded but sparse as well, 

The sparseness is due to the non-consecutive numbering of 
the nodes surrounding the yth node. The optimum compact 
storage (OCS) scheme compacts the matrices by storing only 
the non-zero elements of the matrices. The implementation 
of the OCS scheme requires two additional integer arrays, 
say JA and NAME, The NAME array identifies the nodal points 
which contribute to each nodal equation. The JA array acts 
as a pointer to indicate where the nodal aquaticn starts in 
NAHE. Consider the following simple 2x2 system with nodes 


as indicated 


The NAHE array starts with node 1 and identifies the nodes 
which contribute to node 1 (i.e., 5, 4, and 2). The NAME 


array would, then, give the nodes contrihuting to node 2 


and so forth, so that 
& « $678915 


123 #91 
<NAME> = <1542 : 254162 
i 3 
<JA> = <] 5 V1 oases? 
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3652 °... } 9856 


and 


The algorithm to assemble the element matrices into a compact 


storage vector is straightforward and represents a significant 
savings in computer storage [8]. The system matrices are 
stored as Gc vector rether than a two-dimensional array. -For 
example, the value which would be stored in position (1,5) 


of the nxn array is stored in position 2 of the system vector. 
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VI. NUMERICAL SOLUTION 


This section contains e brief description of possibie 
solution techniques in addition to the solution technique 
chosen. Computer subroutines necessary to implement the 


technique are also described. 


A. SELECTION OF METHOD 

The numerical solution of the system of implicit ordi- 
nary differential equations, equations (64) and (78), may 
be accomplished by any of a number of different technicues 
such as Houbolt's method, Crank-Nicolson's method, Gear's 
method, or implicit Gear's method. It was not the objective 
of this analysis to determine which of the numerical solu- 
tion schemes is the most efficient. Each method has its 
advantages and disadvantages. The Crank-Nicolson sethod 
is a single-step, implicit equation solver and, therefore, 
does not require storage of previous time solutions. When 
analyzing neutronic problems, the system of equations which 
arises is commonly very stiff (i.e., a rapid change in flux 
over a short period of time). The Crank-Nicolson methed 
has, in a past work [8], demonstrated difficulty in tracking 
these stiff systems, Gear's method was specifically 
developed for stiff systems and can handle the probles very 
well, However, Gear's methcd is a multi-step, predictor- 
corrector method requiring storage of previous tise solutions. 


In addition to this disadvantage, Gear's method requires the 


transformation cf the devziuped implicit 0.D.E.'s into an 

explicit system of 0.D.E.'s. After this transformation is 
done, the s:otem matrices are no longer sparse or banded, 

thus eliminating the use of the optimum compacting scheme. 
In an effort to overcome these difficulties, Gear's method | 
was mou.fied, Ref. [9], to treat the implicit system of 

equations &$ well as to allow use of the optimum compacting 


scheme. A previous work, Ref. [8], has shown that the | 


implicit Gear's method is particularly attractive in solv- | 
ing the type problem daveloped in this analysis. Therefore, 

the implicit Gear‘s nethod is used for the solution of the 

system of 9.D.E.'s arising in this analysis. 

No attempt wi?] be made hore to give the mathematics in- | 
yelyed in developing the implicit Gear's method. Reference : 
(9] may be consulted if details are desired. A listing of | 
the computer program developed will be given in the Computer 
Program soction. In order to utflize the implicit Gear's 
method, sevara? user suppiied subroutines must be developed: 
1) DIFFUN, 2) JACMAT, and 3) HUITSL. 

B. USER SUPPLIED SUBROUTINES TO IMPLEMENT 
THE IMPLICIT GEAR'S METHOD 
1, DIFFUN 
Subrout'ne DIFFUR evaluates equations (64) and (78) 
for a givan time and for given values of ¥, }, tc andt,. 


Since at each nodal point, i, there is a solution for the 


flux and for the temperature, the solution was set equal to 


DYI and DYii, respectively. In addition to having flux and 


temperature at each nodal point, there are also three differ- 
ent regions in the domain which have different governing 
equations. An-integer array, ITYPE, was developed to indi- 
cate for each nodal point whether it was: 0) a fuel node 
not in an interface element, 1) a fuel node in an interface 
element, 2) a cladding node or, 3) a coolant node, Using 
ITYPE, she computer program is directed to a different sec- 
tion depending upon the type of node being considered, After 
all the nodes have: baen considered, boundary conditions are 
established by changing DYI and DYII for the appropriate 
boundary nodes, Since in this analysis there is continuity 
of flux at the interfaces, special considerations must be 
given to these nodes. At the fuel-clad interface, the value 
‘of DYI for the clad node was set to the value of the flux 

at that node minus the value of the flux at the adjacent 
node (1.@., DYI, = - Vauy)e Similarly, at the clad-cool- 
ant interface, the value of DYI for the coclant node was 

set to the value of tne flux at that node ininus the value 

of the flux at the adjacent node. During the solution of 


the problem, DYI is driven toward zero, which in the limit 
forcas ¥, to equal 


result. 


2, dACMAT 


Wy.y ¢ This is the desired continuity 


Subroutine JACHAT, evaluates the Jacobian matrix 
(for Gear's method) at the given time and for the current 
values of the dependent variables. The Jacobian for an equa- 
tion of the type, 


Fly, Ys t) « 0 


(80) 


where ao and 8, are coefficients from Gear's method and 
h is the time step. Using the notation of DIFFUN, let DY! 
and DYII represent equations (64) and (78), respectively. 


The Jacobian matrix may, then, be written as {J is called 
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PW in JACMAT.) 

a , a 

PW = poovl o adY!I a0YI! adYTI 

It is the form of equation (81) which is programmed in | 
JACMAT. As in DIFFUN, the integer array ITYPE is used to | 
indicate the appropriate section of the program to be | 
utilized. The problem boundary conditions must also be 
accounted for in JACMAT. In DIFFUN, the value of DYI or 


DYII was set to zero for constant boundary conditions (1.e., 
zero). This cannot be done in JACMAT since « division by 


zero would occur, For a constant boundary condition at the 


ith node, the value of PW is set to one for the diagonal 
term and zero for all other terms of the ith equation. 
3. NUITSL 
Subroutine NUITSL solves the system of equations 
for the quasi-Newton iterates. In this analysis the system 


1s solved using a successive over-relaxgtion (SOR) method. 
In this work, the optimum amount of over-relaxation was not 


determined, Since no effort was made to find the optimum, 


1t was felt a small over-relaxation would be best. The over- 
relaxation factor of 0.02 was used. For small values of 

this factor, the SOR method approaches the Gauss-Siedel 
iteration technique. 
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VII. PROCEDURE 


In this section, the method utilized to obtain a solu- 
tion is described. The input data necessary to run the 
developed computer pregram will be documented. 

Prior to initiating a transient overpower excursion, the 
steady-state conditions for the fuel cell must be known, 
Since the systemof equations which were developed are not 
specifically designed to obtain a steady-state solution, the 
initial steady-state conditions must be part of the input 
data, The initial temperature distribution was obtained 
from the steady-state conditions given in Ref. [7]. 

The axial temperature distribution for the fuel center- 
line, fuel surface, clad, and coolant have been determined 
for several different fuel life cycles [7]. For this analy- 
sis, the beginning of life cycle for channel 10 was used. 
Although this distribution {is somewhat artificial, it should 
be adequate for this analysis. It is the trends of the 
results which are considered important. The distributicn 
within the fuel radially is taken to vary as the square of 
the radial distance; then 


\ nr? re 
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Within the cladding and the coolant, the initial radial tem- 
peratura distribution is assumed to be constant, 
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The initial flux distribution is assumed to be radially 
constant, a flat flux assumption. In the axial direction 
the flux is assumed to vary as the shape of the sine func- 
tion. The maximum flux, the flux at the axial center, is 
an input parameter. For this analysis, the maximum initial 
flux was taken to be 104 neutron/cm“sec. 

To obtain a steady state flux distribution, the value 
of fission cross section for the fuel is varied, A trial- 
and-error method is used until a critical fission cross sec- 
tion, Zep , which gives a steady flux is obtained, 

Once the steady-state conditions have been determined, 
the excess reactivity may be inserted, This starts the 


transient overpower excursion, 


A. INPUT DATA 

The first data card contains: the order of Gauss quadra- 
ture, the number af radial elements in the fuel, the number 
of axial elements, number of nodal points in the radial 
direction, and the height of the fuel rod. The next cards, 
one for each radial nodal point, contain the nodal radial 
distances. The next cards contain the fuel centerline, fuel 
surface, clad and coolant temperaturas. There is one card 
for each axial node, The next card contains the maximum 
flux. The next four ca.ds contain the physical parameters 
listed in Table I, The last input cards contain the time, 


end time, estimated initial time step, minimum time step, 


and maximum time step. A sample data deck is shown in figure 
7. 
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TABLE I. Physical Parameters 


Fuel Diffusion Coefficient (DCF) 
Doppler constant (B) 

Energy release per fission (E) 
Fraction of delayed neutrons (BETA) 


Fraction of delayed neutrons for 
the ith group (BETAI) 


Decay constant for ith: delayed neutron 
group. (DCLAMI 


Initial flux for delayed neutrons 


Average number of neutrons released 
per fission (ANU) 


Neutron velocity (VEL) 

Fual absorption cross section (SIGAF) 

Critical fuel fission cross section 
(SIGFF) 


Blanket absorption cross section (SIGAB) 
Blanket fission cross section (SIGFB) 
Step reactivity input (RHOA) 

Ramp reactivity input (RHOB) 

Fuel density (DENF) 

Clad diffusion coefficient (DCC) 

Clad specific heat (CPC) 

Clad density (OENC) 

Clad thermal conductivity (TKC) 

Clad absorption cross section (SIGAC) 
Coolant diffusion coefficient (DCCO) 
Coolant absorption cross section (SIGACO) 
Coolant flow velocity (VCO) 

Surface heat transfer soeffictent(HSURF) 


0.93 Lom) 
0.006 
7,652x10"'* (cal/#issions] 
0.0064 


0.0064 


0.0784 a 
1xi0!° [neutron/em*sec] 


2.44 

4.8 x10° [cm/sec] 
0,088 Com] 
0,0586875 [om™"] 


0.0 Lem™!y 

0.0 Len"! 
Variable 

Variable 

10.9 [gn/cn?] 

1.1 [om] 

0.12 [cat/gm °C] 
8.0 [gm/em?] 
0.0526 [cal/em sec °C] 
0.0015 Lom™!J 

1.55 [em] 

0.00004 [en™!} 
396.0 [cm/sec] 

0.7 [cal/em*sec °C] 
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VITI. RESULTS 


When using the finite element method, one of the first 


considerations must be given to the convergence of the 


method. To determine convergence, the results for a given 


i 
f 
, 


point are compared for different finite element discretiza- 
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tions. As shown in figure 8, the results are comparable 


mgr 


but no definite claim of cenvergence can be made. However, 
for this work, it was felt that these results were adequate. 
It was not the object of this analysis to arrive at the 
“final” result; it was the trends and methods that were of 
interest. Since the 66-element mesh appears to give a fair 
approximation of the results, the 66-element mesh was used 
as the discretized domain. 

The next item of consideration was the determination of 
a neutronic steady-state condition. This proved to be a 
very time consuming task. The fission cross section for the 
fuel was varied by a trial-and-error method in an attempt 
to find the critical cross section which would give a steady 
state. As may be seen in figure 9, a change in cross sec- 
tion of less than one percent significantly affected the 
state of the problem. It was felt that the critical value 
was between the values of 0.05875 and 0.058625. Time did 
not permit investigation for critical value; therefore, 


it was assumed that the value for the critical fission cross 


section was half way between the values (i.e., 2°°= 0,.0586875), 
f 
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Even if this value is in error, which {it most probably jis, 
the net effect would only be a smal] decraase or increase 
in the proposed reactivity insertion. The reactivity inser- 
tion would, then, be an approximation of the actual reactiv- 
ity of the problem. 
The Tirst test problem considered was a step increase 

in reactivity of approximately ten dollars. For the uranium 
dicxide fuel, one dollar of reactivity was taken to be 
0.0064. Figure 10 shows the time history of the flux at 

the center of the fuel rod. Figure 1] gives the correspond- 
tng temperature profile. The temperatures werc taken at the 
hottest point of each region at the axial center (i.e., the 
fuel centerline, clad inside surface, and coolant inside 
surface). As seen on the fuel temperature time history, the 
fuel rapidly reaches the fuel melting point. The model 
developed does not take into consideration melting of the 
fuel. This melting would tend to decrease the effect of the 
transient. The problem was allowed to continue despite this 
inconsistency in the mathematical model. A short time after 
fuel melting, the inside surface of the clad reaches its 
Melting point. The temperature in the coolant experienced 
What is felt to be a numerical phenomenon. The coolant tenm- 
perature decreased prior to the small rise at the end of the 
transient. Intuitively, this decrease does not seem to be 
realistic, A similar occurrence was observed while conduct- 
ing sample tests on the developed computer program. Refer- 


ence[20] reported the same phenomenon. It is felt that this 
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phenomenon ist quirk of the finite element method. The 
radial and axial distributions of the neutron flux are pre- 
sented in figures 12 and 13 for time equal to 7.39 x 107° 
Seconds. The distributions are, basically, as anticipated. 
The neutron flux peaks stightly before the axial canter. 

It was expected to peak at the axiai center. ‘The radial 


and axial temperature distributions for the same time are 


presented in ficures 14 and 15. As with the flux, the tem- 


perature profiles were, basically, as expected. The fuel 


temperature peaks slightly below the expected location, 
most likely in. response to the peak in axial fiuxes. The 
coolant unexpectedly drops near the sutlet of the fuel rod, 
The finite element method characteristically has some prob- 
tems on the boundaries ef the domain; this may account for 
the drop in coolant temperature. 

The temperature af the fuel and cladding do nat appear 
to be as closely coupled as anticipisted. As seen in figure 
14, a significant increase in fuel temperature has resulted 
ih a relatively small clad temperature increase. As noted 
in figure 11, there appears to be a time lag in temperature 
response for each region which may account for part of the 
apparent temperaturo disagreement, It is felt that the 
temperatures should be more closely related, which indicates 
» higher gap heat transfer coefficient should be utilized. 
Since values of the gap heat transfer coefficient were 
ussumed, it is set unreasonable to belfeve the values used 


are too low. 
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Figure 12. Radial Flux 
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Figure 14. Radial Temperature Profile 
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insertions. Other reactivity inputs may be investigated by 


Time did not permit investigation of other reactivity 
/¢ students in the future. 


This work daes not represent a solution to the very con- 
Slicated nuclear reactor problem. It does represent an 
application of a numerical technique which is relatively new 
to nuclear applications. Methods for implementing the finite 
element method have been discussed, and a computer code has 


been developed for the simplistic model considered, 


IX. RECOMMENDATIONS 


For the model developed, perhaps the most important item 
to pursue is the critical fission cross secticn. A better 
determination of this value is necessary so shat the reactiv- 
ity insertion is more accurately known. Different test cases 
for the prompt critical and prompt subcritical reactor cculd 


then be conducted, 


In further developing the model, more consideration 


should be given the gap heat transfer coefficient. As noted 
in the results, the value used appears to be too small. 
Sample problems for different gap heat transfer coefficients 
would give a better indication of the values to use, 

Melting of the fuel during the transient would prosshly 
be the next major improvement on the model. With relatively 
few changes, the model could be adapted to allow melting 
elenent by element. This, too, would be an approximation 
but, still, an improvement to the model, Perhaps at the 
same time, a simplified model to take into consideration the 
fuel restructuring could de implemented, 

Another improvement would be to consider reactivity 
feedbacks in addition to the Doppler feedback. Sodium void- 
ing and fuel rod expansion are two of the more important 
feedback effects to consider, 

On the numerical side, probably the snost important thing 
to do would be to run the computer program on the 


"Hecompiler", which optimizes the program. However, on 
several runs using the H-compiler, erroneous results were 
Obtained. With sufficient time, this could be corracted 

to allow use of the H-compiler. The use of the H-compiler 
results ina savings in computer time. The present program 
runs on a “Gecompiler" and takes excessive amounts of com- 
puter time (two to four hours per run). 

In addition to this, the optimum over-relaxation factor 
in the implicit Gear's method could be determined by trial- 
and-errcs. 

Implementation of these recommendations should enhance 


the analysis and lead to a more efficient computer code. 


APPENDIX A 
DEVELOPMENT QF TRANSFORMATIONS 


The Jacobian matrix [J] may be written (8) for two 


dimensions as 
N i 
DL Ny vets Nae? 
is] i=] 


(Jj = (Al) 


N N 
Ds Westy Da Neti 


i=] {1 


For a simple 2x2 matrix [A] the inverse is 
a ) 
c d 


(A: ao |: "| 
etlAs =¢ a 


ee ee a ee 


eS a ee 


From matrix algebra 


N N N N 
det[J] > Ser »S Ne ned > Sueur 2 He 24 (A3) 
is] {=1 i=] i=] 


The derivatives of the shape functions may be found fron 


equations (38). 


Nie * q(l4n) (2€4n) Nn T( 146) (2nt6) 

Noe = -E(1ea) Noy = y(tee) 

Kye * glen) (2E=n) Ny, * P18) (2n-8) 

Nee "7 $(1-n*) Ng, =~ nli-é) (A4) 
Neg glen) (2é#n) Hg lle8) (2nt6) 

Neg * ~8(Ton) Neon t  gtleee) 

Nye * glen) (26-n) Kyo = PC148) (2n-8) 

Neg = g(l-n”) Ney = 7n(148) 


I? one now considers an arbitrary element 
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ee ee ee 
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J 
I 
’ 
t 
Fr 


ba} 


1 #86F2 F4 i 


with the midside nodes exactly at the midpoint (not a neces- 
sary criteria for the FEM) and substitutes into equation (A3), 


the result i3 


(zg-Z)(rg-ry) ae 
——— 2 oe 


det[J] = m7 


(A5) 


Substituting into equation (A2) and using (A5) will yield 


4, * °C [F(z5-2,)] = 2/rg-r, »  (A6) 
Jo = (A7) 
a (a8) 


and Joo 2 °C [g(rg-r,)] = 2/232, ; (A9) 


The inverse of the Jacobian matrix now becomes 


tay"! - (A10) 


For integration along a line, the transformation used 


{1s 
az = det[J']dn (All) 
in this case 
N 
det(u'] = DOM, 24 (A12) 
is} 


Again considering the arbitrary element and substituting 
(A&) into (A12) will result in 


det{aty = Sob = Lt (a13) 


eth lll. 


| APPENDIX B | 
| REDUCTION OF SECOND ORDER TERM ! 


7 | The second order term of the governing field equations 
may be reduced to first order by integration by parts. 


Consider, for example, 
Sf ws NL 12 (ro 3 ty +i (0 2 oy) rdrdz (B1) 


which may be expanded as 


2 x 2 
f 3 a9 3 ay a .. . 30 
Sf [N,r 2 +r ue + N,D ae + Np + tHe 4) drdz 
r Zz 


(82) 


a Integrating just the second order terms by parts will yield 
: 2 
a a 
Sf nyo =$ rardz = f [Nyrd SEI dz 
rz z i 
3 20 er 

~ ff setOn, + 1H, 5+ TD x] drdz (83) 
rz 


and 


Sf x0 4 rdrdz dl (N,D zl, rdr S Yi crea DN, a rérdz {B4) 
rz 


Substituting the results of (B3) and (84) into equation (82) 


will give 


feo 3 sel dz + ne Fl rer J i Eat 3 a rdrdz (85) 


8. 
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APPENDIX C 
LIST OF RELATIONS FOR MATERIAL THERMAL PROPERTIES 


FUEL (U0,) 
1. Specific Heat, Ref. £19] 
Cop * [18.45+2.431%10757-2,272x107°121/270.07 
[cal/gm °C] 
T= °C 
2, Thermal Conductivity, Ref. [5] 


4§.1 


Tke = (1-2.8(lepz9) IxEyqese + 4,79x107)? 


79)x0,239 
[cal/cm sec °C] 

Pry ~ percent theoretical density 

T= °K 


CLAD (Stainless Steel) 
Properties are assumed to be temperature independent, 
and average values from Ref. [$] were used for the 
clad properties. 
COOLANT (Liquid Sodium) 
1, Specific heat, Ref. [5] 
0.34574-0.79226x107*1+0,34086x107-T* 
Ccal/gm °C] 
T + °F 
2. fansity, Ref. [5] 
Poo * [89.566-7.9504x107 51-0, 2872x107 872 
+ 0,06035x107°T?; x 0.01601 [om/cm?] 
T = °F 


3. Theraal Conductivity, Ref. [5] 
Teg * (54,306-1.878x10"“T+2.0914x10" 
x 6.234x107° [cal/cm sec *¢} 


T- °F 
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